\documentclass[12pt, a4paper]{article}
\usepackage{geometry}
\geometry{
	a4paper,
	left=12.7 mm,
	right=12.7 mm,
	top=12.7 mm,
	bottom=12.7 mm,
}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{physics}
\graphicspath{ {./} }
\usepackage{ctex}

\newcommand{\bvec}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\formula}[1]{\text{式} \ref{#1} }
\begin{document}
	\section{牛顿第二定律的Verlet算法}
	
	\subsection{引入：牛顿第二定律}
	接下来我们关注一类特殊的二阶ODE，牛顿第二定律。牛顿第二定律可能是最重要的物理定律之一：
	\begin{equation}
		\bvec F=m \bvec a = m\dv[2]{\bvec r}{t}
	\end{equation}
	为了简明起见，我们先考虑一维的牛顿第二定律：
	\begin{equation}
		x=x(t),
		\qquad
		\left \{
		\begin{aligned}
			\dv[2]{x}{t} &= \frac{F(x)}{m} \\
			x(t=t_0) & = x_0 \\
			x'(t=t_0) & = v_{0} \\
			t \ge & t_0 \\
		\end{aligned}
		\right.
	\end{equation}
	其中
	$m$ 是粒子质量，
	$x$是粒子位置，
	$\dv[2]{x}{t}  = a$是粒子加速度，
	$F$是粒子所受力等。
	作为一个二阶ODE，我们不仅要知道粒子的初始位置，还要知道粒子的初始速度(初始的导数条件)。
	
	据我所知，牛顿第二定律之所以特殊，是因其内含特殊的辛结构（通俗地说，能量守恒），
	而之前的高阶ODE迭代法无法保持它的辛结构（会导致较明显的能量误差积累），最终无可避免地导致发散。
	\textsl{不要问我为什么，我也想知道。}
	
	\subsection{Verlet算法}
	
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.5\linewidth]{verletiter}
		\caption{Verlet 算法示意图}
		\label{fig:verletiter}
	\end{figure}
	
	以下先介绍一种简单直观的推导，我们先使用数值微分的方法离散化二阶导 $\dv[2]{x}{t}$：
	\begin{equation}
		\dv[2]{x}{t} |_{t=t_k}  = \frac{x^{(k+1)} - 2x^{(k)} + x^{(k-1)}}{(\Delta t )^2 } \\
	\end{equation}
	那么牛顿第二定律化为
	\begin{equation} 
		\begin{aligned}
			& \dv[2]{x}{t}  |_{t=t_k}  = \frac{F(x^{(k)})}{m} \\
			\Rightarrow & \frac{x^{(k+1)}- 2x^{(k)} + x^{(k-1)}}{(\Delta t )^2 } = \frac{F(x^{(k)})}{m} \\
			\Rightarrow & x^{(k+1)} = 2x^{(k)} - x^{(k-1)}  +\frac{F(x^{(k)})}{m} (\Delta t)^2 \\
		\end{aligned}
		\qquad k=0,1,2,3,...
	\end{equation}
	或者改一下写法
	\begin{equation} 
		x^{(k)} = 2x^{(k-1)} - x^{(k-2)}  +\frac{F(x^{(k-1)})}{m} (\Delta t)^2 \qquad k = 1,2,3,...\\
	\end{equation}
	我们发现，计算 $x^{(k)} $时得知道前两时刻的值$x^{(k-1)} $与$x^{(k-2)} $，
	这就导致一个问题：当我们取$k=0$计算$x^{(1)}$时，我们需要用到$x^{(-1)}$，即$t_0 - \Delta t$时刻粒子的$x$坐标，而这个值好像是不知道的，并且有违反题设之嫌。
	好在，我们可以用初始速度估算$x^{(-1)}$：
	\begin{equation}
		\begin{aligned}
			x^{(-1)} &= x(t_0 - \Delta t)   \\
			& = x(t_0) - \dv{x}{t}|_{t=t_0} \Delta t + \frac{1}{2} \dv[2]{x}{t}|_{t=t_0}  (\Delta t)^2 \qquad \text{取Taylor展开}\\
			& = x_0 - v_{0} \Delta t + \frac{F(x_0)}{2m} (\Delta t)^2 \\
		\end{aligned}
	\end{equation}
	不显式使用$x^{(-1)}$的方法参考附件。总结一下：
	\begin{equation} \label{eq_secder_x}
		\boxed
		{
			\begin{aligned}
				x^{(k+1)} &= 2x^{(k)} - x^{(k-1)}  +\frac{F(x^{(k)})}{m} (\Delta t)^2 \\
				x^{(-1)} &= x_0 - v_{0} \Delta t + \frac{F(x_0)}{2m} (\Delta t)^2 \\
			\end{aligned}
			\qquad k =0,1,2,3,...
		}
	\end{equation}
	尽管思路上很简单，却误打误撞地得到了正确的迭代格式。
	这个方法还称Verlet算法
	\footnote{关于Verlet算法，https://www.zhihu.com/question/22531466 有一个很精彩的讨论}，
	Verlet在他的原始论文
	\footnote
	{
		Loup Verlet. 
		Computer "Experiments" on Classical Fluids. 
		I. Thermodynamical Properties of Lennard-Jones Molecules. 
		Phys. Rev. 159, 98 – Published 5 July 1967
	}
	中使用了这种格式的迭代。
	
	Verlet算法常被用于求解多粒子体系的运动，比较著名的运用是\textsl{分子动力学}。
	分子动力学假定原子的运动遵循经典力学规律，并通过计算各个原子的运动来预测物质的性质，
	因此就需要大量使用Verlet算法。 \textsl{分子动力学挺好玩的，大家可以试一试。}
	
	\subsection{速度Verlet算法}
	
	\begin{figure}[h]
		\centering
		\includegraphics[width=0.5\linewidth]{frogleap}
		\caption
		{
			速度Verlet算法迭代示意图。
			1. 计算$v^{(k+1/2)}$，即半步长的速度；
			2. 根据$v^{(k+1/2)}$等计算$x^{(k+1)}$，即下一时刻的位置;
			3. 根据$x^{(k+1)}$等计算$v^{(k+1)}$，即下一时刻的速度。
		}
		\label{fig:frogleap}
	\end{figure}
	
	\footnote{参考：Peter Young, The leapfrog method and other “symplectic” algorithms for integrating Newton’s laws of motion}
	或许，更优雅的Verlet算法是速度Verlet算法。
	速度Verlet算法原则上等价于上述算法（将在附录中证明），
	但避免使用了$x^{(-1)}$、能同时导出粒子在各个时刻的速度并改善了数值精度等。
	速度Verlet算法使用了定义在半步长的速度：
	\begin{equation} 
		\boxed
		{
			\begin{aligned}
				v^{(k+1/2)} &= v^{(k)} + \frac{F(x^{(k)})}{m} \Delta t / 2 \\
				x^{(k+1)} &= x^{(k)} +v^{(k+1/2)}  \Delta t  \\
				v^{(k+1)} &= v^{(k+1/2)} + \frac{F(x^{(k+1)})}{m} \Delta t / 2 \\
			\end{aligned}
			\qquad k =0,1,2,3,...
		}
	\end{equation}
	如果不喜欢 $v^{(k+1/2)}$，那么我们可以将其合并入其余两项，得到：
	\begin{equation} \label{eq_vv}
		\begin{aligned}
			x^{(k+1)} &= x^{(k)} +v^{(k)} \Delta t  + \frac{1}{2} \frac{F(x^{(k)})}{m} (\Delta t)^2\\
			v^{(k+1)} &= v^{(k)} + \frac{1}{2} (\frac{F(x^{(k)})}{m} + \frac{F(x^{(k+1)})}{m}) \Delta t\\
		\end{aligned}
		\qquad k =0,1,2,3,...
	\end{equation}
	
	
	\subsection{三维Verlet算法}
	在直角坐标系下，三维下的牛顿第二定律由三个形式相似的标量方程组成：
	\begin{equation} \label{eq_n2}
		\bvec F = m \bvec a 
		\Rightarrow
		 \dv[2]{\bvec r}{t} = \frac{\bvec F}{m}
		 \Rightarrow
		 \left \{
		 \begin{aligned}
		 	\dv[2]{x}{t} &= \frac{F_x(x,y,z)}{m} \\
		 	\dv[2]{y}{t} &= \frac{F_y(x,y,z)}{m} \\
		 	\dv[2]{z}{t} &= \frac{F_z(x,y,z)}{m} \\
		 \end{aligned}
		 \right.
	\end{equation}
	其中，$m$是粒子质量，$\bvec r$是粒子坐标（位矢），$ \bvec a = \dv[2]{\bvec r}{t} $是粒子加速度，$ \bvec F$是粒子所受外力；
	此处以$F_x, F_y, F_z$表示分量，而非求偏导。
	
	由于力的形式$\bvec F= \bvec F(x,y,z)$十分多样，这下我们很难解析求解了。
	即使是最好的数学家，对于多体问题的解析解（三体足矣），也爱莫能助。
	不过，我们可以对每一个标量方程分别运用Verlet算法：
	\begin{equation}
		\left \{
		\begin{aligned}
			x^{(k+1)} &= 2x^{(k)} - x^{(k-1)}  +\frac{F_x(x^{(k)},y^{(k)},z^{(k)})}{m} (\Delta t)^2 \\
			y^{(k+1)} &= 2y^{(k)} - y^{(k-1)}  +\frac{F_y(x^{(k)},y^{(k)},z^{(k)})}{m} (\Delta t)^2 \\
			z^{(k+1)} &= 2z^{(k)} - z^{(k-1)}  +\frac{F_z(x^{(k)},y^{(k)},z^{(k)})}{m} (\Delta t)^2 \\
		\end{aligned}
		\right.
	\end{equation}
	压缩为向量形式：
	\begin{equation}
		\bvec r^{(k+1)} = 2 \bvec r^{(k)} - \bvec r^{(k-1)}  + \frac{\bvec F(\bvec r^{(k)})}{m}  (\Delta t)^2  \qquad k = 0,1,2,3,...
	\end{equation}
	估算$\bvec r^{(-1)}$的方法也是类似的。
	
	\newpage

	\section{附件}
	
	\subsection{不显式使用$y^{(-1)}$}
	
	根据导数的定义，$y^{(0)}$处导数$y'_0$可以由$y^{(1)}$与$y^{(-1)}$的函数值估计，由此
	\begin{equation}
		y'_0 = \frac{y^{(1)} - y^{(-1)}}{2 \Delta t} \qquad y^{(-1)} = y^{(1)} - 2y'_0 \Delta t
	\end{equation}
	将$y^{(-1)}$的这个表达式代入$k=0$时的递推关系式 \formula{eq_secder_x} ，得到
	\begin{equation} \label{eq_secder_nonx-1}
		\begin{aligned}
			y^{(1)} &= 2y^{(0)} - (y^{(1)} + 2y'_0 \Delta t)  + f(y^{(0)}) (\Delta t)^2 \\
			\Rightarrow y^{(1)} & = y^{(0)} + y'_0 \Delta t  + 1/2 f(y^{(0)}) (\Delta t)^2
		\end{aligned}
	\end{equation}
	也就是说，如果你不想显式使用$y^{(-1)}$，那么你就得单独计算 $y^{(1)}$：
	计算$y^{(1)}$时使用 \formula{eq_secder_nonx-1}，
	而计算$y^{(2)}, y^{(3)}, y^{(4)}, ...$时再使用原先的递推关系  \formula{eq_secder_x}。
	
	这两种做法实则没有区别，回顾使用$y^{(-1)}$的方法：
	$$
	\left \{
	\begin{aligned}
		y^{(1)} &= 2y^{(0)} - y^{(-1)}  + f(y^{(0)}) (\Delta t)^2 \\
		y^{(-1)} & = y^{(0)} - y'_0 \Delta t + 1/2 f(y^{(0)}) (\Delta t)^2
	\end{aligned}
	\right.
	$$
	得到
	$$
	\begin{aligned}
		y^{(1)} &= 2y^{(0)} - (y^{(0)} - y'_0 \Delta t + 1/2 f(y^{(0)}) (\Delta t)^2)  + f(y^{(0)}) (\Delta t)^2 \\
		& = y^{(0)} + y'_0 \Delta t  + 1/2 f(y^{(0)}) (\Delta t)^2\\
	\end{aligned}
	$$
	与  \formula{eq_secder_nonx-1} 是完全一样的。
	

	
	\subsection{证明Verlet算法的两种形式相同}
	我们将证明\formula{eq_secder_x}等效于\formula{eq_vv}。
	复读一遍\formula{eq_vv}是：
	\begin{equation} \label{eq_ap_3_1}
		x^{(k+1)} = x^{(k)} +v^{(k)} \Delta t  + \frac{1}{2} \frac{F(x^{(k)})}{m} (\Delta t)^2\\
	\end{equation}
	\begin{equation} \label{eq_ap_3_2}
		v^{(k+1)} = v^{(k)} + \frac{1}{2} (\frac{F(x^{(k)})}{m} + \frac{F(x^{(k+1)})}{m}) \Delta t\\
	\end{equation}
	因此，
	\begin{equation} \label{eq_ap_3_3}
		x^{(k)} = x^{(k-1)} +v^{(k-1)} \Delta t  + \frac{1}{2} \frac{F(x^{(k-1)})}{m} (\Delta t)^2
	\end{equation}
	\begin{equation} \label{eq_ap_3_4}
		v^{(k)} = v^{(k-1)} + \frac{1}{2} (\frac{F(x^{(k-1)})}{m} + \frac{F(x^{(k)})}{m}) \Delta t
	\end{equation}
	\formula{eq_ap_3_4}代入\formula{eq_ap_3_1}
	\begin{equation}
		\begin{aligned}
			x^{(k+1)} &= x^{(k)} +(v^{(k-1)} + \frac{1}{2} (\frac{F(x^{(k-1)})}{m} + \frac{F(x^{(k)})}{m}) \Delta t)\Delta t  + \frac{1}{2} \frac{F(x^{(k)})}{m} (\Delta t)^2\\
			& = x^{(k)}+ v^{(k-1)}\Delta t +  \frac{1}{2} \frac{F(x^{(k-1)})}{m}  (\Delta t)^2 +  \frac{F(x^{(k)})}{m} (\Delta t)^2 \\
		\end{aligned}
	\end{equation}
	由 \formula{eq_ap_3_3}
	\begin{equation}
		v^{(k-1)} \Delta t  = x^{(k)} - x^{(k-1)} - \frac{1}{2} \frac{F(x^{(k-1)})}{m} (\Delta t)^2
	\end{equation}
	因此
	\begin{equation}
		\begin{aligned}
			x^{(k+1)} &= x^{(k)}+ x^{(k)} - x^{(k-1)} - \frac{1}{2} \frac{F(x^{(k-1)})}{m} (\Delta t)^2 +  \frac{1}{2} \frac{F(x^{(k-1)})}{m}  (\Delta t)^2 +  \frac{F(x^{(k)})}{m} (\Delta t)^2 \\
			& = 2 x^{(k)}  - x^{(k-1)} +  \frac{F(x^{(k)})}{m} (\Delta t)^2 \\
		\end{aligned}
	\end{equation}
	与 \formula{eq_secder_x} 相同，得证。
	
	\subsection{证明Verlet算法与常规高阶ODE迭代不同}
	将常规高阶ODE迭代改写为只含$x$的：
	\begin{equation}
		\text{常规高阶ODE迭代:} \qquad x^{(k+1)} = 2 x^{(k)}  - x^{(k-1)} +  \frac{F(x^{(k-1)})}{m} (\Delta t)^2 \\
	\end{equation}
	与Verlet算法 \formula{eq_secder_x} 对比
	\begin{equation}
		\text{Verlet算法:} \qquad x^{(k+1)} = 2 x^{(k)}  - x^{(k-1)} +  \frac{F(x^{(k)})}{m} (\Delta t)^2 \\
	\end{equation}
	二者结构上相似，但是Verlet使用了$(k)$步的力，而非$(k-1)$步的力。
	
	
\end{document}
